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Abstract 

Development of an incompressible Navier-Stokes solution procedure was performed for 
the analysis of a liquid rocket engine pump components and for the mechanical heart assist 
devices. The solution procedure for the propulsion systems is applicable to incompress- 
ible Navier-Stokes flows in a steadily rotating frame of reference for any general complex 
configurations. The computer codes were tested on different complex configurations such 
as liquid rocket engine inducer and impellers. As a spin-off technology from the turbop- 
ump component simulations, the flow analysis for an axial heart pump was conducted. 
The baseline Left Ventricular Assist Device (LVAD) design was improved by adding an 
inducer geometry by adapting from the liquid rocket engine pump. The time- accurate 
mode of the incompressible Navier-Stokes code was validated with flapping foil experi- 
ment by using different domain decomposition methods. In the flapping foil experiment, 
two upstream NACA 0025 foils perform high-frequency synchronized motion and generate 
unsteady flow conditions for a downstream larger stationary foil. Fairly good agreement 
was obtained between unsteady experimental data and numerical results from two differ- 
ent moving boundary procedures. Incompressible Navier-Stokes code (INS3D) has been 
extended for heat transfer applications. The temperature equation was written for both 
forced and natural convection phenomena. Flow in a square duct case was used for the 
validation of the code in both natural and forced convection. 


1. Introduction 

Computational fluid dynamics has been developed to a level where it has become an 
indispensable part of aerospace research and design, advanced component systems in the 
next generation of transport vehicles, such as advanced liquid rocket engine fuel systems, 
marine propulsion systems, and high-lift configurations, are likely to require more efficient 
and simpler designs with lower manufacturing costs. Fast and reliable incompressible flow 
simulation capabilities are crucial to developing these new designs. For example, in order to 
understand the source of the noise generated by a high-lift wing-flap system, the unsteady 
flowfield around these configurations must be resolved accurately. 


The cover illustration shows the particle traces colored by relative velocity magnitude in the liquid 
rocket engine inducer, the pressure contours for the Space Shuttle Main Engine HFFTP Impeller, and 
the pressure contours for the old and new designs of the NASA-JSC/Baylor College of Medicine Left 
Ventricular Assist Device (LVAD). 
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Second example of a viscous incompressible flowfield is the rotor-stator interaction in tur- 
bomachinery. A flange-to-flange rocket engine fuel-pump simulation must account for the 
interaction between the rotating and non-rotating components: the flow straighteners, the 
impeller, the volute and diffusers. Such incompressible flow simulation capability can also 
be beneficial to other applications, such as the development of artificial heart devices. A 
Ventricular Assist Device (VAD) developed by NASA Johnson Space Center (JSC) and 
Baylor College of Medicine (BCOM) has a design similar to a rocket engine fuel pump in 
that it also consists of a flow straightener, an impeller, and a diffuser. Technology devel- 
oped for aerospace applications can be utilized for the benefit of human health. Accurate 
and detailed knowledge of the flowfield obtained by steady and unsteady incompressible 
flow calculations can be greatly beneficial to designers in their effort to reduce the cost 
and improve the reliability of these devices. 

Analysis of biofluid flow problems axe challenging and the potential impact of an 
advanced computational tool on medical research for improving health care is enormous. 
For example, the blood flow through mechanical hearts, heart assist devices, and heart 
valves is unsteady and involves moving walls making experimental study very difficult. 
Computational analysis offers an alternative to experimental methods which can produce 
flow field data in great detail. Medical researchers can study this extensive data to obtain 
a better understanding of the flow physics. Computational analysis can also be used 
to optimize the design of mechanical devices at a significantly lower cost and time than 
required by an empirical approach. 

In addition to the geometric complexities, a variety of flow phenomena are encoun- 
tered in biofluids. These include turbulent boundary layer separation, wakes, transition, 
tip vortex resolution, three-dimensional effects, and Reynolds number effects. In order 
to increase the role of Computational Fluid Dynamics (CFD) in the design process, the 
CFD analysis tools must be evaluated and validated so that designers gain confidence in 
their use. The incompressible flow solver, INS3D, has been applied to various viscous flow 
problems and extensively validated. 2-13 This report presents a continuation of this valida- 
tion effort for the turbulent flow in a liquid rocket engine pump components in a steadily 
rotating frame of reference, for time-accurate calculations with moving boundaries, and 
for incompressible flows with heat transfer. In the report, it is also presented how com- 
putational flow simulation capability developed for liquid rocket engine pump component 
analysis has been applied to the Left Ventricular Assist Device being developed jointly by 
NASA JSC and Baylor College of Medicine. 

2. Algorithm 

The algorithm is based on the pseudocompressibility method as developed by Chorin. 14 
The pseudocompressibility algorithm introduces a time-derivative of the pressure term into 
the continuity equation; the elliptic-parabolic type partial differential equations are trans- 
formed into the hyperbolic-parabolic type. The original version of the INS3D code 2-3 with 
pseudocompressibility approach utilized the Beam- Warming 1 5 approximate factorization 
algorithm and central differencing of the convective terms. Since the convective terms of 
the resulting equations are hyperbolic, upwind differencing can be applied to these terms. 
The current versions of the INS3D code use flux-difference splitting based on the method 
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of Roe. 16 Chakravarthy 1 ' outlines a class of high-accuracy flux-differencing schemes for 
the compressible flow equations. Following the third-order upwind differencing, a fifth- 
order-accurate, upwind-biased stencil was derived by Rai. 18 The third- and fifth-order 
upwind differencing used here is an implementation of these schemes for the incompress- 
ible Navier-Stokes equations. The upwind differencing leads to a more diagonal dominant 
system than does central differencing and does not require a user-specified artificial dis- 
sipation. The viscous flux derivatives are computed by using central differencing. In the 
steady-state formulation, the time derivatives are differenced using the Euler backward 
formula. The equations are solved iteratively in pseudo-time until the solution converges 
to a steady state. In the time-accurate formulation, the equations are iterated to conver- 
gence in pseudo-time for each physical time step until the divergence of the velocity field 
has been reduced below a specified tolerance value. 

The matrix equation is solved iteratively by using a nonfactored Gauss-Seidel type 
line-relaxation scheme, 19 which maintains stability and allows a large pseudo-time step 
to be taken. Details of the numerical method can be found in Refs. 4-5. The present 
calculations use the one-equation turbulence model developed by Baldwin and Barth. 20 
In this model, the transport equation for the turbulent Reynolds number is derived from 
a simplified form of the standard k — e model equations. The model is relatively easy 
to implement because there is no need to define an algebraic length scale. The transport 
equation is solved by using the same Gauss-Seidel type line-relaxation scheme as the mean- 
flow equations. 


3. Liquid Rocket Engine Pump Components 

Until recently, the high performance pump design process was not significantly differ- 
ent from that of three decades ago. During this time, a vast amount of experimental and 
operational experience has demonstrated that there are many important features of pump 
flows that are not accounted for in the current semi-empirical design process. Pumps being 
designed today are no more technologically advanced than those designed for the Space 
Shuttle Main Engine (SSME). During that same time span, huge strides have been made 
in computers, numerical algorithms, and physical modeling. One major accomplishment 
of this work is to extend CFD technology by validating an advanced CFD code for pump 
flows and demonstrating their value to pump designers. 

The validation effort for pump components has focused on inducers and impellers. 9-10 
The computed results from a numerical simulation of the flow through a rocket pump in- 
ducer were compared with experimental measurements. The design of an advanced impeller 
was verified by simulating the flow through the baseline and the optimized geometries. The 
effects of the vaneless space in the advanced impeller concept were investigated by using 
the incompressible Navier-Stokes flow simulation capability. The comparisons reported in 
Ref. 9-10 showed that the computations compared well with experiments. The resulting 
computational procedure with the one-equation Baldwin-Barth turbulence model was ap- 
plied to the flow through the SSME High Pressure Fuel Turbo-Pump (HPFTP) impeller. 
The SSME-HPFTP impeller has 6 full blades, 6 long partial blades, and 12 short partial 
blades. The impeller wheel speed is 6,300 rpm, and the Reynolds number for this calcu- 
lation is 181,000 per inch. Figure 1 shows the nondimensional pressure contours on the 


3 


hub surface of the impeller. Since the pressure rise between the impeller inlet and the 
exit is quite large, the pressure variations between the suction and the pressure sides of 
the blades are not noticable from this figure. This calculation includes the vaneless space 
between the impeller trailing edge and the diffuser vane leading edge. Figure 2 shows 
blade-to-blade velocity and flow angle distributions at two stations downstream from the 
SSME-HPFTP impeller exit. Blade-to-blade velocity distribution illustrates the impeller 
exit flow distortion. Solid lines represent the experimental data and dashed lines represent 
computed results. The jet-wake pattern, which produces an unsteady load on the diffuser 
vanes, was captured at both meridional locations. The numerical results compare reason- 
ably well with the experimental data. The pressure contours in Figure 1 also indicates this 
jet-wake pattern. 


4. Left Ventricular Assist Device 

In 1989, NASA Johnson Space Center (JSC) began a joint project with the DeBakey 
Heart Center of the Baylor College of Medicine (BCOM) in Houston to develop a new 
implantable LVAD prototype system. This LVAD is based on a fast rotating axial pump 
requiring a minimum number of moving parts. To make it implantable, the device has 
been made as small as possible, requiring a very high rotational speed. The computational 
procedure described above for pump components has been used to provide the designers 
with a view of the complicated fluid dynamic processes inside this device. Due to the 
nature of the device, this detailed information cannot be obtained experimentally. The 
detailed computational look at the fluid flow is very important to the designers; high levels 
of turbulence can damage the red blood cells and regions of recirculating flow can lead to 
blood clots. Thus the ability to predict these phenomena can greatly help the designers. 

The flow through the baseline design of the LVAD impeller was numerically simulated 
by solving the incompressible Navier-Stokes equations in a steady rotating frame of refer- 
ence. Zonal multiblock grids were used in this component analysis. The geometry and the 
computational grid of the LVAD baseline impeller are shown in Figure 3. The domain is 
divided into five zones with dimensions of 127 x 39 x 33, 127 x 39 x 33, 59 x 21 x 7, 47 x 
21 x 5, and 59 x 21 x 7, respectively. Zone 1 is the region between the suction side of the 
partial blade and the pressure side of the full blade; the region between the pressure side 
of the partial blade and the suction side of the full blade is filled by zone 2; and zones 3 
through 5 allow tip-leakage effects to be included in the computational study and occupy 
the regions between the impeller blade tip and the casing. At the zonal interfaces, grid 
points were matched one-to-one. For all zones, an H-H type grid topology was used. An 
H-type surface grid was generated for each surface using an elliptic grid generator. The 
interior region of the three-dimensional grid was filled using an algebraic grid generator 
coupled with an elliptic smoother. Periodic boundary conditions were used at the end 
points in the rotational direction. The design flow of this impeller is 5 liters per minute 
and the design speed is 12,600 revolutions per minute (rpm). The problem was nondimen- 
sionalized by the tube diameter (0.488 inches) and the impeller tip velocity. The solution 
was considered converged when the maximum residual had dropped at least five orders 
of magnitude. Computer time required per grid point per iteration was about 1.5 x 10 -4 
seconds. The total computer time required for these calculations was about 6 to 8 single 
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processor Cray-C90 hours. A parametric study was performed to optimize the impeller 
blade shape and the tip clearance. Initially, three different impeller blade designs with a 
tip clearance of 0.009 inches were analyzed. Figure 4 shows the axial velocity contours 
at the exit plane of the impeller for various blade curvatures. Design III shows massive 
seperation near the suction side of the partial blade trailing edge because of the high cur- 
vature blade shape. Design II, which has less blade curvature than design I, shows flow 
features very similar to design I. Design I was analyzed .with two tip clearances; the tip 
clearance of 0.0045 inches shows better hydrodynamic performance in terms of efficiency 
and head coefficient than with a tip clearance of 0.009 inches. Using the design I with a 
tip clearance of 0.0045 inches as the baseline impeller design, ideas from rocket propulsion 
and medical science were combined to develop a new implantable LVAD. In collaboration 
with the NASA-JSC engineering team and BCOM researchers, a new design consisting 
of the baseline impeller plus an inducer was investigated. The hub and blade surfaces of 
the baseline impeller and the new impeller are colored by nondimensionalized pressure in 
Figure 5. The pressure is nondimensionalized by pV 2 , where p is the density and V is the 
impeller tip velocity. The pressure gradient across the blades due to the action of centrifu- 
gal force, and the pressure rise from inflow to outflow are shown. The table in Figure 5 
shows the clinical results obtained by BCOM. Further clinical hemolysis testing and ani- 
mal implantation by BCOM are currently underway using the new inducer-impeller design. 
The hemolysis index reported in Figure 5 shows the amount of hemoglobin generated by 
the pump in grams per 100 liters. Hemoglobin release can be due to regions of high shear 
and separated flows. The new design shows a remarkable improvement in performance over 
the baseline design. The inducer provides a sufficient pressure rise to the flow in order to 
prevent the cavitation on the impeller blades. Figure 6 shows the particle traces through 
the new impeller. The traces are colored by the relative total-velocity magnitude. The 
particles were released near the inducer leading edge, the hub, the inducer blade pressure 
side, and the tip regions. The swirling motion of the particles indicates a secondary flow 
region between the partial and the full blades. The particles released near the pressure 
side of the blade indicate a radial velocity component inside the blade boundary layer. The 
particles tend to flow from the hub to the tip of the blade. The particles near the inducer 
leading edge and full blade trailing edge indicate the presence of a back flow. 

In the next step of the design process, this computational analysis tool will be used 
to improve the diffuser geometry and the bearing design. This unique insight into the 
internal fluid structures will lead to improved components for an improved heart assist 
device. In July 1991, the Institute of Medicine estimated that approximately 25,000 to 
60,000 patients per year in North America could benefit from an efficient left ventricular 
assist device. Thus improved designs made possible because of the current work could 
have a far reaching impact on human health. 


5. Heat Transfer 

Heat transfer in viscous incompressible flow is of great interest in many engineering 
applications. For example, in a liquid rocket engine, the liquid fuel and oxidizer are used 
as coolant in various components, and in the LVAD, the electric motor releases heat to the 
blood flow. In order to predict the effect of the heat transfer, a simplified energy equation 
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is added to current incompressible Navier-Stokes formulation. The fluid is assumed to be 
incompressible with constant physical properties except the bouyancy effect due to density 
variations. Using the Boussinesq approximation for bouyancy force, the incompressible 
Navier-Stokes equations with temperature equation in nondimensional form can be stated 
in the following way: 


du{ 

dxi 

dui 

~dt 

dT 

~di 


+ 
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du. i u 
dxj 
dTuj 
dxj 


L _ dp 
dxi 

= c<v 2 t 


+ C v \? 2 u — e g CbT 


( 1 ) 


where u,- is the corresponding velocity components, p the pressure, T the temperature, t 
the time, and aq the cartesian coordinates. C v , Cb, and C t are nondimensional coefficients 
which are defined as for forced (or mixed) convection, 
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and for natural convection, 


C v = Pr , C b = RaPr , C t = 1 
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Here, u 0 , is the reference velocity, l a the reference lenght, ( T\ — To), is the reference 
temperature difference, u the kinematic viscosity, a. the thermal diffusivity, and /? the 
coeficient of expansion. For natural convection, the reference velocity u Q is replaced by 
a// 0 in the nondimensionalization procedure. 

In order to validate the current heat-transfer formulation with Boussinesq approxima- 
tion, fluid flow with heat transfer in a square cavity is considered. First, forced convection 
problem was solved by setting Gr = 0. In this case the momentum equations and the 
temperature equation are decoupled. A 83 x 83 mesh was used for the Reynolds number 
400. The computed velocity magnitde and temperature contours for this driven-cavity flow 
is ploted in Figure 7. Figure 7 also shows the boundary conditions used for this computa- 
tion. A large temperature gradient is observed near the hot and cold vertical walls. Forced 
convection results show very good agreement with previous study reported in reference 21. 
For the Rayleigh number 10,000, free convection results are presented in Figure 8. In this 
case, the lid remains stationary, and an oposite direction vortex to the one in the forced 
convection case is formed due to heat transfer. The temperature distribution shows very 
good agreement with computed results obtained by velocity-vorticity formulation. 21 Figure 
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9 plots u velocity component along the vertical midplane and v velocity component along 
the horizontal midplane. Along the midplanes, the peak values of the velocity components 
reported in reference 21 were also plotted. Very good agreement is obtained between the 
two. More validation cases for heat transfer will be reported in future. 


6. Time-Accurate Calculations 

The validation effort for pump flows was extended to unsteady interaction problems by 
simulating the flapping foil experiment performed by Massachusetts Institute of Technology 
(MIT). The purpose of the MIT flapping foil experiment (FFX) is to provide detailed 
unsteady measurements for CFD code validations for the marine propulsion systems. The 
FFX was designed as a two-dimensional representation of the interaction between the 
propeller-blade and wake flows. The stationary foil represents a propeller blade embedded 
in a wake flow generated by upstream pitching foils. The flappers in the FFX generate 
periodic wake fluctuations which impose an unsteady condition on the stationary foil. 

First, steady-state calculations were carried out for overlapped and patched grids 
in which effects of the grid density and of the order of differencing were investigated. 
Numerical results showed good agreement with the experimental data. Then, the resulting 
third-order upwind differencing and Baldwin-Barth one equation turbulence model were 
applied to unsteady flapping foil calculations. Fairly good agreement was obtained between 
unsteady experimental data and numerical results from two different moving boundary 
procedures. Appendix B includes the results from these computations. 


7. Conclusions 

An efficient and robust solution procedure for 3-D pump component analyses and its 
spin-off application to an LVAD impeller analysis has been presented. The technique solves 
the viscous incompressible Navier-Stokes equations with source terms in a steadily rotating 
reference frame. The method of pseudocompressibility with higher-order accurate upwind 
differencing and a Gauss-Seidel line relaxation scheme are utilized. The flow through 
SSME-HPFTP impeller has been successfully simulated. Numerical results which utilize 
the one-equation Baldwin-Barth turbulence model compare fairly well with experimental 
data. This validated solution procedure has been applied to the NASA JSC/BCOM LVAD 
impeller as a design analysis tool. Parametric studies were performed to analyze the 
performance of the LVAD impeller. Substantial improvements were observed when inducer 
blades were included upstream of the impeller main blades. 

Steady and time-accurate calculations were performed for the MIT flapping foil ex- 
periment. A detailed study for the steady-state calculations was performed to investigate 
the effects of grid topology, grid resolution, and the order of accuracy. A nearly grid 
independent solution was obtained for the steady FFX simulation by using third-order 
flux-difference splitting for the convective terms. Unsteady computations were performed 
using two different grid topologies. Good agreement with experimental mean pressure 
coefficient were obtained for both grid topologies. 
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The temperature equation is coupled with momentum equations in order to include 
heat transfer effect in incompressible flow calculations. Using the Boussinesq approxima- 
tion for bouyancy force, driven cavity problem was studied for forced and mixed convection 
case. Very good agreement is obtained between the computed results and the results re- 
ported in reference 21. 
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Figure 2 : Velocity and flow angle distributions downstream of 
SSME-HPFTP impeller exit plane (impeller exit radius : 5.5 in.) 
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Figure 3 : Geometry and computational grid for the old design 
of the LVAD impeller 
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Figure 4 : Axial velocity contours at the exit plane of various LVAD impeller design 
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Figure 5 : Hub and blade surfaces of the old and new designs of the LVAD 
impeller are colored by pressure. The table shows clinical results from BCOM 
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Figure 6 : Particle traces inside of the new design of the LVAD 
impeller arc colored by velocity magnitude 
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Figure 7 : Velocity vectors and temperature contours (forced 
convection. Re = 400 ) 
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Figure 8 : Velocity vectors and temperature contours (natural 
convection, Ra = 10,000 ) 



Figure 9 : Profiles of u- and v- velocity components along 
vertical and horizontal mid-planes (Ra = 10,000) 
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Abstract 

Steady and unsteady flows for propulsion sys- 
tems are efficiently simulated by solving the incom- 
pressible Navier-Stokes equations. The solution method 
is based on the pseudocompressibility approach and 
uses an implicit-upwind differencing scheme together 
with the Gauss-Seidel line relaxation method. Current 
computations use one-equation Baldwin-Barth turbu- 
lence model which is derived from a simplified form 
of the standard k — e model equations. The resulting 
computer code is applied to the flow analysis inside 
an advanced rocket pump impeller in steadily rotat- 
ing reference frames. Numerical results are compared 
with experimental measurements. The effects of exit 
and shroud cavities with the leakage flow are inves- 
tigated. Time-accurate incompressible Navier-Stokes 
formulation with the overlapped grid scheme capabil- 
ity was evaluated by using MIT flapping foil exper- 
iment. The effects of grid density and of the order 
of differencing were investigated. The resulting pro- 
cedure were applied to unsteady flapping foil calcula- 
tions. Two upstream NACA 0025 foils perform high- 
frequency synchronized motion and generate unsteady 
flow conditions to the downstream larger stationary 
foil. Numerical results obtained from steady flow com- 
putations were compared against experimental data. 

Introduction 

Since the space launch systems in the near future 
are likely to rely on liquid rocket engines, increasing 
the efficiency and reliability of the engine components 
is an important task. One of the major problems in 
the liquid rocket engine is to understand the fluid dy- 
namics of fuel and oxidizer flows. Understanding the 
flow in the turbopump through numerical simulation 
will be of significant value toward finding a better de- 
sign that is simpler yet more efficient and robust with 
less manufacturing cost. Until recently, the pump de- 
sign process was not significantly different from that of 
decades ago. The current semi-empirical turbomachi- 
nary design process does not account for the three- 
dimensional (3-D) viscous phenomena in the pump 
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flows. Some of these 3-D viscous phenomena include 
wakes; the boundary layers in the hub, the shroud and 
the blades; junction flows; and tip clearence flows. In 
order to meet the challenge of improving propulsion 
devices, NASA Marshall Space Flight Center (MSFC) 
has established a consortium involving universities, in- 
dustries, and NASA. 1-2 

Even though computational fluid dynamics (CFD) 
applications in turbines have been reported widely in 
the literature 3-5 , the applications in the pump area 
are quite limited. The objective of this paper is to 
present, evaluate, and validate a computational pro- 
cedure that solves incompressible Navier-Stokes equa- 
tions for the pump components. For CFD to have an 
impact in the design procedure of the pump, a robust, 
efficient, and accurate scheme is required. In addi- 
tion, the algorithm needs to be extensively validated 
for the flow-through pump components so that pump 
designers have the confidence to use it. The present 
work is focused on steady-state component analyses to 
validate the algoritm with a one-equation turbulence 
model and to demonstrate the code capability. For the 
pump components, such as an inducer and a radial im- 
peller, the steady flow assumption is valid without the 
diffuser and inlet guide vanes. The progress in un- 
steady pump flows will be reported in the future. An 
effort similar to the present work can be seen in the 
activities of the other members of the MSFC Pump 
Stage Technology (PST) team. 6-8 Their formulations 
are based on a pressure based method and the cur- 
rent formulation is based on a pseudocompressibility 
method. 

The numerical solution of the incompressible 
Navier-Stokes equations requires special attention in 
order to satisfy the divergence-free constraint on the 
velocity field because the incompressible formulation 
does not yield the pressure field explicitly from the 
equation of state or through the continuity equation. 
One way to avoid the numerical difficulty originated 
by the elliptic nature of the problem is to use a pseu- 
docompressibility method. With the pseudocompress- 
ibility method, the elliptic-parabolic type equations 
are transformed into hyperbolic-parabolic type equa- 
tions. Well established solution algorithms developed 
for compressible flows can be utilized to solve the re- 
sulting equations. Steger and Kutler 9 employed an 
alternating-direction implicit scheme into Chorin’s 10 


pseudocompressibility method. This formulation was 
extended to three-dimensional generalized coordinates 
by Kwak and Chang. 11,12 Recently, a three-dimensional 
incompressible Navier-Stokes solver 13 (INS3D-UP) that 
uses upwind differencing and the Gauss-Seidel line re- 
laxation scheme was developed in order to have a ro- 
bust and fast-converging scheme. The flow over sin- 
gle and multi-element airfoils w'as simulated efficiently 
by using this algorithm. 14 A time-accurate formula- 
tion of the algorithm is implemented for incompress- 
ible flows through artificial-heart devices with moving 
boundaries. 15 The present study is a continuation of 
the validation effort for the turbulent flow in rotating 
machinery in a steadily rotating frame of reference. 
Preliminary results from a benchmark case for the ma- 
rine propulsors are also included in this paper. 

In the following sections, the governing equations 
and the method of solution are summarized, and the 
computed results obtained from the current approach 
are presented. 

Algorithm 

The pseudocompressibility algorithm introduces 
a time derivative of the pressure term into the conti- 
nuity equation. The resulting incompressible Navier- 
Stokes equations can be written in a generalized curvi- 
linear coordinate system 33 follows 

where Q and the convective flux vectors E, F , G are 
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Here 7, /?, p, u, v, and w denote the Jacobian of trans- 
formation, the pseudocompressibility coefficient, pres- 
sure, and cartesian velocity components, respectively. 
The contravariant velocity components (7, V, and W 


are defined as 
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The viscous fluxes, E v , F u , and G „, are given by 
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where Re is the Reynolds number. 

When the equations are solved in steadily rotat- 
ing reference frames, the centrifugal force and the Cori- 
olis force terms are added to the equation of motion as 
source terms. If the relative reference frame is moving 
around the x— axis, the source term S is given by 

r 0 

s - ® 

Q(Qy -f 2 w) 

. — 2v) . 

where Q is the rotational speed. Both the viscous and 
source terms are treated implicitly in the numerical 
formulation. Relative velocity components are written 
in terms of absolute velocity components u 0 , u„, and 
w a as 

v = u a 
v = v a + firr 
w = ui 0 — Qy 

In the steady-state formulation, the time derivatives 
are differenced using by the Euler backward formula. 
The equations are solved iteratively in pseudo-time 
until the solution converges to a steady state. Cen- 
tral differencing is used to compute the viscous flux 
derivatives and third-order upwind differencing is em- 
ployed to compute the convective flux derivatives. 
Chakravarthy 16 outlines a class of high-accuracy flux- 
differencing schemes for the compressible flow equa- 
tions. Following Chakravarthy’s third-order scheme, 
a fifth-order-accurate, upwind-biased stencil was de- 
rived by Rai. 3 The upwind differencing used here is an 
implementation of those efforts for the incompressible 
Navier-Stokes equations. For the computations pre- 
sented in this paper, third-order flux- 
differencing scheme is used. 

An implicit, delta-law form approximation to equa- 
tion (1) after linearization in time and the use of ap- 
proximate Jacobians of the flux differences result in a 
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seven-block diagonal matrix equation written as 

B6Qi-i,j,k + ASQij'k + C6Q{+ ij-.jt -f DSQij-.i't 
+E6Qij + ik + FSQij't^i + GSQij't+i = R.H.S. 

M V) 

where 6Q = Q n+l - Q n and A,B,C t D,E,F, and G 
are 4x4 block diagonals. 

The Gauss-Seidel line relaxation scheme, which 
was successfully employed by MacCormack 17 , is used 
to solve this matrix equation. In equation (2), the 
right-hand-side term is computed and stored for the 
entire domain. In the present study, the line-relaxation 
procedure is composed of three stages; each stage in- 
volves a block-tridiagonal inversion in one direction. 
Fig 1 shows the block tridiagonal inversion lines and 
the Gauss-Siedel sweep planes for the three directions. 

In the first stage, 6Q is solved line by line in one direc- 
tion. Before the block-tridiagonal equation is solved, 
ofT-tridiagonal terms are multiplied by the current value 
of SQ and then shifted over to the right-hand-side of 
the equation. In other words, equation (2) is solved 
by performing a block-tridiagonal inversion in the 
£— direction and Gauss-Seidel sweeps in the q— and 
(^—directions. The second stage is to solve the block- 
tridiagonal terms in the rj— direction, and to perform 
backward and forward sweeps in the (— and directions. 
The same procedure is repeated in the third stage by 
inverting the block-tridiagonal matrix in the ( -direction, 
and treating the off-diagonal terms for the £— and 
^—directions in Gauss-Seidel fashion. After the first 
sweep is completed for the entire domain, a backward 
sweep is started in the opposite direction. One for- 
ward sweep and one backward sweep for each compu- 
tational direction are sufficient for most problems, but 
the number of sweeps can be increased. 

Implicit boundary conditions are used at all the 
boundaries except the zonal interface boundaries; zonal 
interface boundaries are updated quasi-implicitly, a 
method very suitable for the line-relaxation scheme. 
The change in the dependent variables for one time- 
step is passed during line-relaxation sweeping from 
one zone to another. Nonslip boundary conditions are 
imposed at the solid stationary wall. The pressure 
boundary condition is specified such that the pres- 
sure gradient normal to the wall is zero. At the inflow 
and outflow boundaries, characteristic boundary con- 
ditions are employed implicitly. It is assumed that the 
effect of viscous terms at the boundaries is negligible. 
Also, the characteristic equations are approximated in 
one-dimensional space. The primitive variables that 
are needed at the boundaries are the pressure p and 
the velocity components u, v, w. The number of pos- 
itive and negative eigenvalues of the Jacobian matrix 
of the convective flux determines how many variables 
should be specified at the boundaries. If the flow is in 
the positive £— direction, then there will be three char- 
acteristic waves traveling downstream and one char- 
acteric wave traveling upstream. The exit boundary 


receives the information about the three variables via 
the characteristics traveling from the interior of the do- 
main. Hence, only one variable is specified at the out- 
flow boundary. However, the inflow boundary receives 
only one characteristic traveling from the interior re- 
gion, and, therefore, three variables are specified at 
the inflow boundary. For the calculations presneted in 
this report, u, v and w velocities were specified at the 
inflow and static pressure was specified at the outflow. 
Details of the numerical method are given in reference 
13. 

The present calculations use the one-equation tur- 
bulence model developed by Baldwin and Barth. 18 The 
transport equation for the turbulent Reynolds num- 
ber is derived from a simplified form of the standard 
k — e model equations. The model is relatively easy 
to implement because there is no need to define an al- 
gebraic length scale. The formulation and code issues 
can be found in reference 18. The transport equation 
is also solved by using a Gauss-Seidel type line relax- 
ation scheme. 

Pump Components 

The fiowfield through a turbopump inducer ge- 
ometry as shown in figure 2 was solved as a benchmark 
problem for turbomachinery applications. An inducer 
that provides a sufficient pressure rise to the flow in 
order to prevent the cavitation on impeller blades is 
a crucial element of a rocket engine pump. The de- 
sign flow of the Rocketdyne inducer is 2236 gallons per 
minute (gal/min) and, the design speed is 3600 revo- 
lutions per minute (rpm). In the computational study, 
tip-leakage effects are included with a tip clearance of 
0.008 inches (in). The Reynolds number for this calcu- 
lation was 191,800 per inch. The solution was consid- 
ered converged when the maximum residual dropped 
at least five orders of magnitude. The convergence 
history is shown in figure 3. Computer time required 
per grid point per iteration was about 1.5 x 10~ 4 sec- 
onds (sec). The surfaces of the inducer hub and blades 
are colored by nondimensionalized pressure in figure 4. 
The pressure is nondimensionalized by pV 2 , where p is 
the density and V is the average inflow velocity. The 
pressure gradient across the blades due to the action 
of centrifugal force and the pressure rise from inflow to 
outflow are shown. Figure 5 shows the particle traces 
colored by relative total-velocity magnitude. The par- 
ticles were released near the hub, the blade suction 
side, and the tip regions, and then the traces were 
calculated from the relative velocity field. Near the 
hub, the particles rise up in the boundary layer and 
move from the pressure side to the suction side of the 
blade. The swirling motion of the particles indicates 
a secondary flow region near the hub. The particles 
released near the suction side of the blade indicate a 
radial velocity component inside the blade boundary 


layer because low energy fluid is controlled by cen- 
trifugal force. The particles tend to flow from the hub 
suction side to the tip suction side of the blade. The 
particles near the casing show an opposite trend from 
the ones near the hub. Since the casing has counterro- 
tating wheel speed in the relative frame of reference, 
the particles move from the suction side to the pres- 
sure side. The structure of the internal turbulent flow 
in the present configuration is quite complicated. The 
comparison between numerical results and exerimen- 
tal measurements was presented in reference 19. The 
comparison showed that the solution algorithm does a 
reasonably good job. The existing solution procedure 
can be applied to a similar configuration in off-design 
conditions. Such a numerical study could potentially 
predict cases in which the inducer may suffer from 
massive separation and result in a blocked fuel sup- 
ply. The numerical study can provide the designer a 
safe operating envelope of a particular inducer. This 
is the future research area of this study which can be 
used in the predesign and postdesign engineering tool 
in challenging turbomachinery applications. 

The current procedure was applied in a flow anal- 
ysis inside an advanced pump-impeller geometry to 
verify the design. Comparison between the computed 
results and experimental measurments was presented 
in reference 20. The impeller design flow is 1,205 
gal/min with a design speed of 6,322 rpm. The Reynolds 
Number for this calculation was 181,273 per inch. Fig- 
ure 6 shows the computational grid near the hub region 
of the impeller. The results with vaneless space and 
without the vaneless space show the effect of down- 
stream conditions. Figures 7-a and 7-b show the merid- 
ional velocity distribution at the impeller discharge. A 
relative x-distance is measured from the shroud to the 
hub, where x = 1.0 is the hub. The meridional veloci- 
ties, CM, were integrated along a radial strip for each 
constant x — position and they were nondimensional- 
ized by the wheel tip speed of 249.5 ft/sec. In figure 
7-a, the dashed line represents CM distributions for 
the flow without vaneless space. The solid line denotes 
the CM distributions for the flow with vaneless space. 
The effect of downstream boundary conditions can be 
seen by comparing the solid line with the dashed line 
in figure 7-a. In case where the vaneless space is not 
included, the flow is pumped near the hub and shroud 
regions, and the velocity profile is flattened in the core 
region. The meridional velocity distributions for 5% 
and 10% recirculation from the exit shroud cavity were 
also plotted in figure 7-b. When the shroud cavity has 
leakage to the impeller eye, the velocity peak at the 
impeller exit moves toward to the center of the B2 
width. However, the effect of the shroud leakage has 
minor effects on the solution at r/rtip = 1.0275 (fig- 
ures 7-b and 8). In figure 7-b, the symbols represent 
experimental data 21 , and the lines represent CM dis- 
tributions for the flow with vaneless space. When the 
vaneless space is included, the velocity profile shows 


a peak at approximately 60% of the B2 width. The 
test data shows that the peak is closer to the center 
of the B2 width. The discrepency between the com- 
puted results and experimental data is partially due 
to the recirculation flow in the hub cavity. The leak- 
age at the hub cavity leads to stronger recirculation 
region which shifts the velocity peak to the center of 
B2 width. Since the CFD analysis did not include the 
leakage at the hub cavity, the predicted recirculation 
region in the vaneless space is not as strong as in the 
experimental study. Figure 8 shows blade-to-blade ve- 
locity distributions at the impeller exit. The jet-wake 
pattern was captured at all three axial locations. The 
numerical results compare fairly well with the experi- 
mental data. 

MIT Flapping Foil Experiment 
Part-I: Steady Flow 

The purpose of the MIT flapping foil experiment 
(FFX) is to provide detailed unsteady measurements 
for CFD code validations for marine propulsors. The 
FFX was designed as a two-dimensional representa- 
tion of the interaction between the propeller-blade and 
wake flows. 22 A schematic of the experimental set-up 
is given in figure 9. The stationary foil represents a 
propeller blade embedded in a wake flow generated 
by upstream pitching foils. The stationary foil has 
an 18-inch chord and a 1.18 degree an gle-of- attack. 
The upstream flapping foils are NACA 0025 foils of 
3-inch chord. The flappers perform synchronized si- 
nusoidal motions of 6-degree amplitude at a reduced 
frequency of 3.62. The flappers in the FFX generate 
periodic wake fluctuations which impose an unsteady 
condition on the stationary foil. Velocity and pressure 
measurements at a Reynolds number (based on the 
stationary foil chord and the inflow mean velocity) of 
3.78x10® were taken on and around the stationary foil 
inside the measurement box, which is shown by the 
dashed line in figure 9. Previous computations of the 
FFX include the following: a simplified configuration 
for the FFX was numerically simulated using a poten- 
tial flow formulation; 23 a patched, time- varying, multi- 
block grid system was employed with a pseudocom- 
pressibility formulation by researchers at Mississippi 
State University 24 to simulate the FFX configuration; 
a chimera overlaid-grid scheme with a pressure-based 
method was used in the FFX simulation by Peterson 
and Stern. 25 In this paper, the results from steady flow 
calculations are presented. Since experimentally mea- 
sured data was available for steady flow with station- 
ary flappers with zero degree angle of attack, steady- 
state calculations were carried out to validate the cur- 
rent approach before starting the more time consum- 
ing time-dependent calculations. Unsteady flapping 
foil calculations will be presented in reference 26. 

The grid topology for the FFX is shown in fig- 
ure 10. An H-type grid with the dimension of 253x191 


(grid 1) occupies the water tunnel without considering 
the foils. Three C-type grids are generated for the foils 
and are overlapped with the tunnel grid. Grid 2 is gen- 
erated for the stationary foil with the grid dimension 
of 337x61. Grids 3 and 4 wrap around the flappers 
with the grid dimension of 215x40. The total num- 
ber of grid points for this grid system is 86080. The 
individual grids receive information from each other 
by interpolating the dependent variables. The grids 
around the foils have outer boundaries which overlap 
the interior region of the tunnel grid. The tunnel grid 
has an interior boundary surrounding a hole. A hole 
point is a mesh point which is removed from the solu- 
tion procedure. The immediate neighbor points of the 
hole points are called fringe points and are updated 
from the interpolation procedure. 

Figure 11 compares the measured and calculated 
static pressure coefficient (Cp) on the stationary foil 
surface. The symbols represent the experimental mea- 
surements, and the lines represent the computed re- 
sults. All numerical results show very good agree- 
ment with the measured data. Figures 12 through 14 
show the velocity profiles at several streamwise loca- 
tions on the surfaces and at the wake of stationary 
foil. These profiles plot the magnitude of the compo- 
nent of velocity that is tangential to the local station- 
ary foil surface. Symbols represent the experimental 
measurements and the solid lines represent numerical 
results obtained by using the overlapped grid topol- 
ogy. Figures 12 through 14 indicate that the computed 
results compare very well with the measured data at 
the boundary layer and the wake of the stationary foil. 
The biggest discrepancy between the computed results 
and the measured data is seen in the wake of the foil 
(x/c=1.2). The velocity at the edge of the wake is 
overpredicted with less then a couple of percent error 
range. In order to make sure the steady flow results 
are grid independent, some additional computations 
were carried out. The order of accuracy for convective 
terms in coarse grid calculations was increased from 
third-order to fifth-order flux difference splitting. The 
third-order coarse grid result is shown with the dashed 
line and the fifth-order coarse grid result is shown with 
the chaindashed line in figure 15. The solid line and 
the dotted line with x-marks represent the fine grid 
third-order and fifth-order results, respectively. The 
wake from the fine grid computations clearly show bet- 
ter agreement than the coarse grid does with the mea- 
sured data. The overshoot which occurs at the edge 
of the wake in the third-order results does not occur 
in the fifth-order results. In the finest grid computa- 
tions, the base grid for the stationary foil was refined 
by increasing the number of grid points in the normal 
direction from 61 to 91. The velocity profile from the 
resulting finer grid (total 97K points) with the third- 
order differencing is plotted with the chaindotted line. 
This fine grid result is very similar to the base grid 
(86K grid point) result, indicating that this is close 


to grid independent solution and that 86K point grid 
will provide adequate resolution for the unsteady cal- 
culations. In fact, the difference between the 86K and 
the 97K grid results is less than the oscillations in the 
measured data. Notethat there is a rather large dif- 
ferencing between the measured and computed wake 
edge velocities. Since this edge velocity is shown to be 
grid independent, it is thought that the experimental 
data used an erroneous value for the reference velocity 
for this data. In addition, all computed results have 
the same velocity magnitude at the edge of the wake 
and fine grid results are self consistent. Because of 
these reasons, the result shown with solid line in figure 
15 is considered to be a reasonably grid independent 
solution. At this point, it is believed that the compu- 
tational procedure is ready for the unsteady flapping 
foil calculations which are currently underway. Re- 
sults from the unsteady calculations will be repoted in 
reference 26. 

Conclusions 

An efficient and robust solution procedure for 
3-D pump component analyses has been presented. 
The technique solves the viscous incompressible Navier- 
Stokes equations with source terms in a steadily rotat- 
ing reference frame. The method of pseudocompress- 
ibility with higher-order accurate upwind differencing 
and a Gauss-Seidel line relaxation scheme are utilized. 
The flow through a pump inducer and an impeller have 
been successfully simulated. Numerical results from 
a one-equation Baldwin-Barth turbulence model com- 
pare fairly well with experimental data. 

A detailed study for the steady-state MIT flap- 
ping foil calculations was performed to investigate the 
effects of grid resolution, and the order of accuracy. 
The one-equation Baldwin-Barth turbulence model was 
able to accurately compute the Cp and boundary layer 
velocity distribution on the stationary airfoil. A nearly 
grid independent solution was obtained for the steady 
FFX simulation by using third-order flux-diflerence 
splitting for the convective terms. Overall, the com- 
puted steady-state results showed good agreement with 
the experimentally measured Cp and boundary layer 
profiles on on the stationary foil. 
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Fig. 1. Gauss-Seidel sweep directions. 



Fig. 2. Rocketdyne turbopump inducer configuration. 



Fig. 3. Convergence history for inducer calculations. 



Fig. 4. Surface pressure of the pump inducer. 



Fig. 5. Particle traces colored by relative total velocity 
magnitude. J 
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Fig. 6. Advanced pump impeller computa- 


tional grid on the hub surface. 
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Fig. 7-a. Effect of the vaneless space on the 
impeller exit velocity. 
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Fig. 7-b. Comparison of circulferentially 
averaged meridional velocity at the impeller 
exit. 
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Fig. 8. Comparison of blade-to-blade merid 
ional velocity at the impeller exit. 
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Fig. 9. Schematic of MIT Flapping Foil 
experiment (FFX). 



Fig. 10. Chimera overlapped grid topology 
for the FFX. 



Fig. 11. Steady-state Cp distribution on the 
stationary foil 



Fig. 12. Velocity profiles on upper surface 
of the stationary foil at streamwise stations 
of x/c=0.388, x/c=0.900, and x/c=0.972. 



Fig. 13. Velocity profiles on lower surface of 
the stationary foil at streamwise stations of 
x/c=0.388, x/c=0.759, and x/c=0.900. 
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Fig. 14. Velocity profiles at the wake of 
the stationary foil at streamwise stations of 
x/c=1.05, x/c=1.10, and x/c=1.20. 



Fig. 15. Velocity profiles at the wake of 
the stationary foil at streamwise station of 
xjc = 1.20. 
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Abstract 

A flapping foil experiment conducted at Mas- 
sachusetts Institute of Technology was used as a val- 
idation case to evaluate the current incompressible 
Navier-Stokes approach with overlapped grid schemes. 
Steady-state calculations were carried out for over- 
lapped and patched grids in which effects of the grid 
density and of the order of differencing were investi- 
gated. Numerical results showed good agreement with 
the experimental data. The resulting third-order up- 
wind differencing and Baldwin-Barth one equation tur- 
bulence model were applied to unsteady flapping foil 
calculations. Two upstream NACA 0025 foils perform 
high-frequency synchronized motion and generate un- 
steady flow conditions for a downstream larger sta- 
tionary foil. Fairly good agreement was obtained be- 
tween unsteady experimental data and numerical re- 
sults from two different moving boundary procedures. 

Introduction 

The next generation of transport vehicles with 
advanced component systems (e.g. high-lift configura- 
tions in the advanced subsonic systems, marine propul- 
sion systems, advanced liquid rocket engine fuel sys- 
tem, etc.) in the near future are likely to require more 
efficient and simpler designs with lower manufacturing 
costs. Therefore, developing fast and reliable time- 
accurate incompressible flow simulation capabilities is 
an important task. For example, in order to under- 
stand the source of the noise generated by a high-lift 
wing-flap system, one needs to accurately resolve the 
unsteady air flow around these configurations. An- 
other example of an unsteady incompressible flowfield 
is rotor-stator interaction in turbomachinary applica- 
tions. A flange-to-flange rocket engine fuel-pump sim- 
ulation includes an interaction between the rotating 
and non-rotating components, such as between flow 
straighteners, inducer-impeller, volute and diffusers. 
The unsteady incompressible flow simulations capabil- 
ity can also be beneficial to artificial heart device devel- 
opers. A Ventricular Assist Device (VAD) developed 
by NASA Johnson Space Center and Baylor College 
of Medicine has flow straightener, inducer-impeller, 


and diffuser interaction which is similar to rotor-stator 
interaction described above for the fuel pump. 1 The 
technology developed for aerospace applications can 
be utilized for the benefit of the human health. Accu- 
rate and detailed knowledge of the flowfield obtained 
by unsteady flow calculations can be greatly benefi- 
cial for designers to reduce the cost and to improve 
the reliability of the advanced systems. In addition to 
geometric complexities, the challenges in these numer- 
ical simulations include turbulent boundary layer sep- 
aration, wakes, transition, tip vortex resolution, three- 
dimensional effects, Reynolds number effects, the time- 
dependency and moving boundaries. As an example of 
complex physics in these examples, preliminary exper- 
imental efforts indicated that small amplitude oscilla- 
tions of vibrating ribbons, flaps or fences delayed Lhe 
separation of turbulent boundary layers and increased 
maximum lift generated by these airfoils. 2 A validated 
time-accurate incompressible Navier-Stokes (INS) so- 
lution procedure can be used to investigate the possi- 
ble ways to control separation whenever the mean flow 
is unstable. In order to increase the role of Compu- 
tational Fluid Dynamics (CFD) in the design process, 
the time-accurate algorithm needs to be evaluated and 
validated so that designers have the confidence to use 
it. 

The current algorithm 3- ' 1 for the incompressible 
Navier-Stokes equations is based on pseudocompress- 
ibility method which was introduced by Chorin. 5 The 
incompressible flow solvers developed at NASA Ames 
Research Center (the INS2D and INS3D family of codes) 
have been applied successfully to aerodynamic, turbo- 
machinary, and artificial heart flow simulations: steady 
flow over high-lift aerodynamic configurations has been 
efficiently and accurately simulated; 6-8 the compo- 
nent analysis of the liquid rocket engine-pump /lows 
validated the current approach in a steadily rotat- 
ing frame of references; 9 wingtip vortex flowfield was 
studied using the present approach in steady-state; 10 
A time-accurate artificial heart flow simulation with 
moving boundaries and some validation cases related 
to blood flow have been computed. 11 The objective of 
the current work is to validate the time-accurate for- 
mulation of the code for the flows with moving bound- 
aries. An ideal validation case for this purpose has 


been supplied by the Office of Naval Research (ONR) 
and Massachusetts Institute of Technology (MIT) Un- 
steady Flow Workshop held on March 29-30, 1993. 
ONR/MIT designed a flapping foil experiment (FFX) 
as a two-dimensional representation of the interaction 
between the propeller-blade and wake flows. 12 One 
purpose of the experiment is to provide detailed ex- 
perimental data to be used to evaluate computational 
methods for marine propulsors. At the ONR/MIT 
workshop, the computed results obtained by various 
group of researchers were compared with experimental 
measurements. The flappers in the FFX generate high 
frequency periodic wake fluctuations which impose an 
unsteady loading on the stationary foil. In addition 
to the complexity of the flow physics, the numerical 
simulation of the FFX requires a proper domain de- 
composition and moving boundary procedures. This 
makes the FFX an attractive validation study for the 
current time-accurate formulation. 

In the following section, the algorithm for the 
solution of incompressible Navier-Stokes equations is 
summarized. Next, the MIT FFX and the computa- 
tional models for the FFX are described. Computed 
results from steady and time-accurate calculations us- 
ing two different moving boundary procedures are then 
presented. 


Algorithm 


The present computations are performed utiliz- 
ing the INS2D-UP 3 computer code which solves the in- 
compressible Navier-Stokes equations for both steady- 
state and unsteady flows. The algorithm is based on 
the pseudocompressibility method as developed by 
Chorin. 5 The pseudocompressibility algorithm intro- 
duces a time-derivative of the pressure term into the 
continuity equation; the elliptic-parabolic type par- 
tial differential equations are transformed into the 
hyperbolic-parabolic type. The original version of the 
INS3D code 13-14 with pseudocompressibility approach 
utilized the Beam- Wanning 1 5 approximate factoriza- 
tion algorithm and central differencing of the convec- 
tive terms. Since the convective terms of the resulting 
equations are hyperbolic, upwind differencing can be 
applied to these terms. The current versions of the 
INS2D and INS3D codes use flux-difference splitting 
based on the method of Roe. 16 Chakravarthy 17 out- 
lines a class of high-accuracy flux-differencing schemes 
for the compressible flow equations. Following the 
third-order upwind differencing, a fifth-order-accurate, 
upwind-biased stencil was derived by Rai. 18 The third 
and fifth-order upwind differencing used here is an 
implementation of these schemes for the incompress- 
ible Navier-Stokes equations. The upwind differenc- 
ing leads to a more diagonal dominant system than 
does central differencing and does not require a user- 
specified artificial dissipation. The viscous flux deriva- 
tives are computed by using central differencing. In 
the steady-state formulation, the time derivatives are 


differenced using the Euler backward formula. The 
equations are solved iteratively in pseudo-time until 
the solution converges to a steady state. In the time- 
accurate formulation, the time derivatives in the mo- 
mentum equations are differenced using a second-order, 
three-point, backward-difference formula. 
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where U , and V are contravariant velocity compo- 
nents, and where q and r denote the dependent vari- 
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Here At, A r, n, m , and /? denote physical time step, 
pseudo-time step, physical time level, subiteration time 
level, and pseudocompressibility coefficient, respectively. 
The equations are iterated to convergence in pseudo- 
time for each physical time step until the divergence 
of the velocity field has been reduced below a specified 
tolerance value. 

The matrix equation is solved iteratively by us- 
ing a nonfactored Gauss-Seidel type line-relaxation 
scheme, 19 which maintains stability and allows a large 
pseudo-time step to be taken. Details of the numerical 
method can be found in Refs. 3-4. The present cal- 
culations use the one-equation turbulence model de- 
veloped by Baldwin and Barth. 20 In this model, the 
transport equation for the turbulent Reynolds number 
is derived from a simplified form of the standard k — e 
model equations. The model is relatively easy to im- 
plement because there is no need to define an algebraic 
length scale. The transport equation is solved by using 
the same Gauss-Seidel type line-relaxation scheme as 
the mean-flow equations. 


MIT Flapping Foil Experiment 
and Computational Models 


The purpose of the MIT FFX is to provide de- 
tailed unsteady measurements for CFD code valida- 
tions. A schematic of the experimental set-up is given 
in figure 1. The stationary foil represents a propeller 
blade embedded in a wake flow generated by upstream 
pitching foils. The stationary foil has an 18-inch chord 
and a 1.18 degree angle-of-attack. The upstream flap- 
ping foils are NACA 0025 foils of 3-inch chord. The 
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Figure 1. Schematic of MIT Flapping Foil ex- 
periment (FFX). 

flappers perform synchronized sinusoidal motions of 6- 
degree amplitude at a reduced frequency of 3.62. Pre- 
vious computations of the FFX include the following: 
a simplified configuration for the FFX was numeri- 
cally simulated using a potential flow formulation; 22 
a patched, time-varying, multi-block grid system was 
employed with a pseudocompressibility formulation by’ 
researchers at Mississippi State University 23 to simu- 
late the FFX configuration; a chimera overlaid-grid 
scheme with a pressure-based method was used in the 
FFX simulation by Peterson and Stern. 24 In the present 
study, various grid topologies for the FFX configura- 
tion are studied, and the computed results obtained 
from a time-accurate pseudocompressibility formula- 
tion are compared with the experimental data. 12 

The flappers in the FFX generate periodic wake 
fluctuations which impose an unsteady condition on 
the stationary foil. Velocity and pressure measure- 
ments at a Reynolds number (based on the stationary- 
foil chord and the inflow mean velocity) of 3.78x10® 
were taken on and around the stationary foil inside the 
measurement box, which is shown by the dashed line in 
figure 1. The box measurements were intended to pro- 
vide the upstream, downstream, and outer boundary’ 
conditions for numerical calculations of the stationary 
foil alone. Since the purpose of the current numerical 
study is to investigate the moving boundary capability 
of the present unsteady’ incompressible Navier- Stokes 
solution procedure, the current computational model 
includes the entire domain shown in figure 1. The ex- 
perimental conditions at the inlet and at the exit of 



Figure 2. Multi-block patched grid topology 
for the FFX. 


the water tunnel are specified numerically. In order to 
simulate the entire configuration, domain decomposi- 
tion methods are utilized. Two commonly used do- 


main decomposition methods for structured grids are 
multi-block patched 15 and chimera overlapped 21 grid 
schemes. 

Figure 2 shows a multi-block patched grid topol- 
ogy' applied to the FFX geometry’, which contains four 
H-type grids. The patched grids are pointwise contin- 
uous at the zonal boundaries, and the interfaces have 
two points of overlap. Each grid has the dimension 
of 319x63, resulting in a total number of 80388 grid 
points. Every other grid lines were plotted in all grid 
related figures in this text. The patched grids were 
generated by GRIDGEN 25 with the elliptic grid gen- 
erator option. Grid 1 occupies the region between the 
lower tunnel wall and the lower flapper surface; grid 2 
extends between the lower flapper upper surface and 
the stationary foil pressure side; grid 3 is located be- 
tween the stationary foil suction side and the upper 
flapper bottom surface; and grid 4 extends between 
the upper flapper top surface and the upper tunnel 
wall. The advantage of having this type of multi-block 
patched grid scheme is that the grids remain pointwise 
continuous as the bodies move, and thus there is no in- 
terpolation error at the interface boundaries. However, 
the grid has to be re-generated at each physical time 
step because of the flappers’ motion. For the FFX, 
the interface boundaries between the zones move up 
or down with the flappers’ rotation, and each zone ex- 
periences contraction and expansion during the cyclic 
motion. 



Figure 3-a. Chimera overlapped grid topology 
for the FFX. 


An alternative to the multi-block patched grid 
scheme is the chimera overlapped grid scheme. The 
overlapped grid topology for the FFX is shown in fig- 
ure 3-a. An H-ty’pe grid with the dimension of 253x191 
(grid 1) occupies the water tunnel without considering 
the foils. Three C-ty'pe grids are generated for the 
foils and are overlapped with the tunnel grid. Grid 2 
is generated for the stationary foil with the grid di- 
mension of 337x61. Grids 3 and 4 wrap around the 
flappers with the grid dimension of 215x40; they ro- 
tate with the flappers. The total number of grid points 
for this grid system is 86080. Grid 1 was generated al- 
gebraically, and the C-type grids were generated by 
the HYPGEN 26 hyperbolic grid generator. The ad- 
vantage of the chimera overlapped grid scheme as a 
moving boundary procedure is the simplification of the 
grid generation procedure. For the FFX, the grids are 
generated once, then the flapper grids are rotated rel- 
ative to the tunnel grid. However, additional numer- 


ical boundary conditions and the data management 
for the time-dependent interpolation stencils are intro- 
duced. The overlapped grid regions in the near-field of 
the flapper and in the near-field of the stationary foil 
leading edge are plotted in figure 3-b. The individual 
grids receive information from each other by interpo- 
lating the dependent variables. The grids around the 
foils have outer boundaries which overlap the interior 
region of the tunnel grid. The tunnel grid has an in- 
terior boundary surrounding a hole. A hole point is 
a mesh point which is removed from the solution pro- 
cedure. The immediate neighbor points of the hole 
points are called fringe points and are updated from 
the interpolation procedure. For all computations pre- 
sented in this paper, two-layer of fringe points are used 




Figure 3-b. Overlapped grid regions in the near- 
field of the stationary foil and the flapper 

for interior and outer boundaries. The PEGSUS 21 
software is used to generate all interpolation informa- 
tion for each discrete position of the flapper grids. This 
grid topology is referred as “the overlapped grid” in 
the later text. For the overlapped grid topology of the 
FFX, the interpolation data between the flapper grid 
and the tunnel grid is time-dependent while the inter- 
polation data between the stationary foil grid and the 
tunnel grid remains constant in time. The major dif- 
ference between the patched and the overlapped grid 
approaches is the amount of effort required to gener- 
ate the time-varying grids. Generating time-varying 


patched grid system for the FFX using the elliptic 
grid generator required an order of magnitude more 
work than generating the overlapped grid system and 
the interpolation data base. The overlapped grid sys- 
tem provides the flexibility of choosing grid topologies 
since the grids do not have constraints at boundaries. 
Therefore, the C-type hyperbolic grids can be easily 
used around the foils in which very fine grid resolu- 
tion is required near the boundary layer. The patched 
grid s 3 'stem designed for the FFX required to use the 
H-type grid with the constraints at the boundaries. 
The elliptic grid generator were used for these Ii-grids. 
Obtaining the desired grid density near the boundary 
layer was the most time-consuming part in this pro- 
cedure which had to be repeated at every boundary 
movement. 



Figure 4. Composite grid topology for the FFX. 


The last grid topology studied in this paper for 
the FFX is illustrated in figure 4. Both patched and 
overlapped grid schemes are employed in this grid sys- 
tem, referred as “the composite grid” in the later text. 
Three H-grids (grids 1,3, and 4) are patched around 
the flappers by using the GR.IDGEN elliptic grid gen- 
erator. A C-grid around the stationary foil (grid 2) was 
generated by using the HYPGEN hyperbolic grid gen- 
erator. In fact, grid 2 for the overlapped grid system 
and grid 2 for the composite grid system are identical. 
A hole is cut in grid 1 to accommodate the station- 
ary foil. Grids 1 and 2 communicate with each other 
through the chimera interpolation procedure. The in- 
terpolation data between the stationary foil grid and 
the tunnel grid is not constant in time because the tun- 
nel grid moves in time. Total number of grid points 
for this composite grid system is 77932, and the di- 
mensions of the grids are 255x99, 337x01, 255x63, and 
255x63, respectively. 

Computed Results 

In this section, the steady flow calculations with 
three different grid topologies are presented. The ef- 
fects of grid density, and of the order of the differenc- 
ing of the convective terms are investigated. Then the 
resulting procedure is applied to unsteady flow calcu- 
lations. 

Steady Flow Results 

Since experimentally measured data was avail- 
able for steady flow with stationary flappers with zero 




degree angle of attack, steady-state calculations were 
carried out to validate the current approach before 
starting the more time consuming time-dependent cal- 
culations. The pseudocompressibility coefficient (/?) 
was taken as 10 for all grid topologies in this sec- 
tion. Figure 5 compares the measured and calculated 
static pressure coefficient (Cp) on the stationary foil 
surface. The symbols represent the experimental mea- 
surements. The dashed line represents the patched 
grid results, the solid line represents the overlapped 
grid results, and the chaindotted line represents the 
composite grid results. All results show very good 
agreement with the measured data. The composite 
grid system Cp results are nearly identical to the over- 
lapped grid results. 



Figure 5. Steady-state Cp distribution on the 
stationary foil 


Total velocity magnitude contours from the over- 
lapped grid calculations are shown in figure 6. The 
wakes of the stationary foil and of the flappers are 
clearly seen. The contours from both grids in the over- 
lapped regions match each other quite well. This en- 
sures that the converged solution for the FFX steady- 
state case were obtained. The convergence history for 
this case is shown in figure 7. Converged solutions 
were obtained after 250 iterations which corresponds 
to approximately six minutes of CPU time on a Cray 
C90. Very similar convergence behavior was observed 
for all grid topologies. Even though the results did not 
change significantly after 250 iterations, 600 iterations 
were performed in order to show the behavior of the 



Figure 6. Steady-state total velocity magnitude 
contours. 


convergence. The solid line in figure 7 shows the his- 
tory of the maximum residual of the mean-flow equa- 
tions. The dashed line shows the history of the maxi- 
mum of the divergence of velocity, and the chaindotted 
line shows the history of the RMS of the Baldwin- 
Barth one-equation turbulence model. Figure 7 indi- 
cates that all three of these measures of convergence 
behave in a similar fashion. 



Figure 7. Convergence history for the FFX 
(overlapped grid topology). 

The velocity profiles from the suction surface 
boundary layer of the stationary foil at the streamwise 
station of x/c = 0.612 are plotted in figure 8 for the 
three different grid topologies. These profiles plot the 
magnitude of the component of velocity that is tan- 
gential to the local stationary foil surface. This figure 
shows the effects of grid resolution near the stationary 
foil wall. The dashed and chaindashed lines represent 
the results from the patched grid system, in which each 
zone has an H-type grid generated by an elliptic grid 
generator. The elliptic grid generator does not provide 
exact control of the grid spacing at solid wall bound- 
aries. Even though spacing on the order of 10~ 5 were 
specified as input to the elliptic grid generator, the re- 
sulting wall spacing was typically on the order of 10 -3 . 
Therefore, the grid resolution near the stationary foil 
wall is quite poor for H-grids in this calculation. The 
velocity profile shown with the dashed line in figure 
8 does not compare well with the experimental data. 
In order to improve the grid resolution near the wall 
region, the elliptic grid generator is being forced not 
to move the near wall points away from the wall. This 
result is shown with the chaindashed line. Although 
the velocity profile shown with the chaindashed line is 
improved compared to the one with the dashed line, it 
still does not compare well with the experimental data. 
The overlapped grid computations for the three levels 
of the grid density were carried out. The total number 
of grid points for the coarse overlapped grid system was 
38607, and the dimensions of the grids were 127x96, 



337x35, 215x34, and 215x34, respectively. The dimen- 
sions of the base overlapped grid system, which had 
total 8G0S0 points, were given in the previous section. 



Figure 8. Velocity profiles on upper surface 
of the stationary foil at streamwise station of 
x/c=0.612. 


ally identical with overlapped grid results (see figure 
8). Figures 9 through 11 indicate that the computed 
results compare very well with the measured data at 
the boundary layer and the wake of the stationary foil. 
The biggest discrepancy between the computed results 
and the measured data is seen in the wake of the foil 
(x/c=1.2). The velocity at the edge of the wake is 



Figure 9. Velocity profiles on upper surface 
of the stationary foil at streamwise stations of 
x/c=0.388, x/c=0.900, and x/c=0.972. 


The total number of grid points for the finest over- 
lapped grid system was 97480, and the dimensions of 
the grids were 253x191, 337x91, 215x43, and 215x43, 
respectively. The dotted line in figure 8 represent over- 
lapped grid results from the coarse grid. The grid 
spacing near the wall for the coarse grid is 2.5xl0~ 4 . 
The dotted line shows better agreement with the ex- 
perimental data than patched grid results. The solid 
line represent the result obtained by using the base 
overlapped grid system and the dashed line with the 
x-maxk symbols represent the result obtained by using 
the composite grid topology. The grid spacing near 
the stationary foil wall for these fine grids is G.OxlO -5 . 
Both velocity profiles are virtually identical and com- 
pare fairly well with the measured data. These calcu- 
lations use C-type hyperbolic grid for the stationary 
foil and clearly show the convenience of that approach 
over the Ii-type elliptic grid for airfoil type of config- 
urations. The amount of work in generating H-type 
elliptic grid was increased when better grid resolution 
was required in the boundary layer region. Because of 
the easiness in the grid generation procedure, it was 
decided to use C-type hyperbolic grid around the sta- 
tionary foil. Therefore, further computations in this 
text were focused on the overlapped grid topology and 
the composite grid topology. It should be pointed out 
that the composite grid topology still has patched H- 
grids around the flappers. 

Figures 9 through 11 show the velocity profiles at 
several streamwise locations on the surfaces and at the 
wake of stationary foil. Symbols represent the experi- 
mental measurements and the solid lines represent nu- 
merical results obtained by using the overlapped grid 
topology. The velocity profiles from the composite grid 
topology are not included here because they are virtu- 



Figure 10. Velocity profiles on lower surface 
of the stationary foil at streamwise stations of 
x/c=0.388, x/c— 0.759, and x/c=0.900. 
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Figure 11. Velocity profiles at the wake of the 
stationary foil at streamwise stations of x/c=1.05, 
x/c=1.10, and x/c=1.20. 


overpredicted with less then a couple of percent error 
range. In order to make sure the steady flow results 
are grid independent, some additional computations 






were carried out. The order of accuracy for convective 
terms in coarse grid calculations was increased from 
third-order to fifth-order flux difference splitting. The 
third-order coarse grid result is shown with the dashed 
line and the fifth-order coarse grid result is shown with 
the chaindashed line in figure 12. The solid line and 
the dotted line with x-marks represent the fine grid 
third-order and fifth-order results, respective!}'. The 



Figure 12. Velocity profiles at the wake of the 
stationary foil at streamwise station of x/c = 
1 . 20 . 


wake from the fine grid computations clearly show bet- 
ter agreement than the coarse grid does with the mea- 
sured data. The overshoot which occurs at the edge 
of the wake in the third-order results does not occur 
in the fifth-order results. In the finest grid computa- 
tions, the base grid for the stationary foil was refined 
by increasing the number of grid points in the normal 
direction from 61 to 91. The velocity profile from the 
resulting finer grid (total 97K points) with the third- 
order differencing is plotted with the chaindotted line. 
This fine grid result is very similar to the base grid 
(86K grid point) result, indicating that this is close 
to grid independent solution and that 86K point grid 
will provide adequate resolution for the unsteady cal- 
culations. In fact, the difference between the 86K and 
the 97K grid results is less than the oscillations in the 
measured data. Notethat there is a rather large dif- 
ferencing between the measured and computed wake 
edge velocities. Since this edge velocity is shown to be 
grid independent, it is thought that the experimental 
data used an erroneous value for the reference velocity 
for this data. In addition, all computed results have 
the same velocity magnitude at the edge of the wake 
and fine grid results are self consistent. Because of 
these reasons, the result shown with solid line in figure 
12 is considered to be a reasonably grid independent 
solution. At this point, it is believed that the compu- 
tational procedure is ready for the unsteady flapping 
foil calculations. 


Unsteady Flow Results 

The converged steady solutions were used as the 
initial conditions for the unsteady calculations. The 
motion of flappers was specified as a pitching about the 
mid-chord point described by a = a m sin(zut), where 
t is time, a,„ is six degree and the angular velocity 
is described by w — 2kUoo/c. Here is the refer- 
ence velocity specified as 20.62 feet/ sec and c is the 
stationary foil chord specified as 18 inches. The re- 
duced frequency k is 3.62 based on c/2. In all unsteady 
calculations, one cycle of the flapping foil period con- 
sisted of 192 physical time steps, which corresponds 
to a discrete nondimensional time step of 4. 53x10. _3 
The time-step size was chosen to be small enough so 
that the moving mesh points in the overlapped grid 
system does not move through more than one cell in a 
neighboring grid during one time step. At each phys- 
ical time-step, the maximum of non-dimensional di- 
vergence of velocity was dropped below 10 -2 for all 
zones. This required between 15 to 40 subiterations 
during each time-step. Numerical tests indicated that 
reducing the divergence of velocity lower than this had 
no effect. The pseudo-time step (A r) was taken the 
same value as physical time step (At). The pseudo- 
compressibility coefficient p was set to 10. With these 
conditions, one period of the flapping foil was obtained 
in 3.5 to 4 Cray C90 hours. The periodic solution was 
obtained for both grid topologies after six flapping 
cycles. 



(a) with the overlapped grid topology 




(b) with the composite grid topology 


Figure 13. Unsteady velocity magnitude con- 
tours obtained by using two different grid topolo- 
gies at t/T = 0.25. 

Total velocity magnitude contours at a nondi- 
mensional time of t/T = 0.25 are shown in figure 13 
for both the overlapped grid topology (top figure) and 
the composite grid topology (bottom figure). Here 
T denotes the period of the flappers’ motion. This 
qualitative comparison shows that the results are very 



similar. This unsteady wake is resolved slightly bet- 
ter in the composite grid topology than it is in the 
overlapped grid topology. Since the grid patch bound- 
aries lay in the region where the oscillating wakes are 
located, the finer grid resolution in the patched grid 
boundary region results in the better resolution of the 
wakes. However, the difference between the two re- 
sults is not so recognizable from the contours. The 
quantitative comparison between the two and the ex- 
perimental data is presented in figures 14 through 16. 



Figure 14. Unsteady mean Cp distribution on 
the stationary foil 


The mean Cp distributions on the stationary foil 
surface from the unsteady calculations are plotted in 
figure 14. The symbols represent the experimental 
mean values. The solid lines represent the overlapped 
grid results, and the dashed lines represent the com- 
posite grid results. The computed mean Cp values 
compare very well with the experimental measurements, 
except there is a slight overprediction at the 60% of 
chord location on the pressure side of the foil. The 
mean Cp values from two different grid topologies show 
identical behavior. The time history of Cp values at 
several streamwise locations of the stationary foil are 
compared with experimental data in figure 15. The 
Cp values obtained from both grid topologies show a 
very similar time history. They both compare fairly 
well with the measured data. The biggest discrepan- 
cies are seen on the pressure side of the foil at stream- 
wise location x/c = 0.611 where Cp is overpredicted 
and on the suction side of the foil at streamwise lo- 
cation x/c = 0.972 where Cp is underpredicted. This 
is consistent with the mean Cp distribution in Fig- 
ure 14. The major difference between the two com- 
puted results is that overlapped grid topology predicts 
some higher frequency oscillations that the compos- 
ite grid scheme does not. It should be noted that the 
grid movement and the number of overset grid bound- 
ary interpolations in the composite grid topology is 
considerably smaller compared to the overlapped grid 


topology. In the overlapped grid topology, as the flap- 
per grids rotate they move through relatively coarse 
region in the tunnel grid. This mismatch of the grid 
resolution between overlapped grids can lead to inter- 
polation errors and an inability of grid to resolve the 
flapper wakes. Considering that different hole points 
are being cut at each time step while locally refined 
flow feature moves in time (unsteady wakes from flap- 
pers), obtaining the same degree of accuracy between 
the two grids is not an easy task. In the composite 
grid system, the mesh points in the tunnel grid move 
with the flappers which maintains relatively fine grid 
resolution in the near wakes of the flappers. Improve- 
ments of the accuracy for the overset grids are cur- 
rently being researched. Recently, a defect correction 
approach 27 has been introduced to enhance accuracy 
for overset grids. The FFX simulation with overlapped 
grid topology would be a good test case for the defect 
correction scheme in the future. 
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Figure 15. Time history of Cp at various stream- 
wise locations of the stationary foil. 


Figure 16 illustrates the effect of the improper 
use of the interpolation information when computing 
previous time level data for newly created fringe points. 
The time history of Cp on the suction side of the foil 
at streamwise location x/c = 810 is replotted for the 
composite grid topology. The symbols represent the 





experimental measurements. The procedures for up- 
dating previous time level data at the fringe points are 
different between the two computed results shown by 
the solid line and by the dashed line. As shown in 
equation 1, information at time levels n and n — 1 is 
neccesarry to advance to the n + 1 time-level. When 
a hole point becomes a fringe point due to moving 
boundaries this fringe point does not have any infor- 
mation from the previous time levels. One way to 
obtain the previous-time level data is to interpolate 
the variables from the donor grid (as it is done for 
the current time- level information). The result ob- 
tained by this procedure is plotted by the dashed line 
in figure 16. This shows that very large amplitude 
errors are occuring in the computations. The source 
of these fluctuations can be found in the interpolat- 
ing procedure. The previous time-level data for the 
new fringe points is obtained by using the interpola- 
tion data base at the current time level. However, 
the grid point locations from the previous time-level 
should have been used. Using the current time-level 
data base for these points results in incorrect interpo- 
lation coefficients and incorrect donor points. When 
this error was corrected, the computed Cp value shown 
by the solid line in figure 16 was obtained. When the 
previous time-level data is not available for the newly 
created boundary points, the time integration for these 
points is changed to first-order. The time-differencing 
in the next time step will be second-order backward 
differencing because the previous time-level informa- 
tion is established from the current time step. 



Figure 16 . The effect of updating boundary 
points in the overlaid moving regions of the 
composite grid system. 


Conclusions 


An incompressible flow solver in steady and time- 
accurate formulations has been utilized for the MIT 
flapping foil experiment. A detailed study for the 


steady-state calculations was performed to investigate 
the effects of grid topology, grid resolution, and the 
order of accuracy. The one-equation Baldwin-Barth 
turbulence model was able to accurately compute the 
Cp and boundary layer velocity distribution on the 
stationary airfoil. A nearly grid independent solution 
was obtained for the steady FFX simulation by us- 
ing third-order flux-difference splitting for the convec- 
tive terms. Overall, the computed steady-state results 
showed good agreement with the experimentally mea- 
sured Cp and boundary layer profiles on on the sta- 
tionary foil. 

Unsteady computations were performed using two 
different grid topologies. Good agreement with ex- 
perimental mean pressure coefficient were obtained for 
both grid topologies. The time history of Cp on the 
stationary foil at several streamwise locations was pre- 
sented for the two grid topologies. Initial attempts 
of the computations on both types of grid topologies 
resulted in high-frequency errors of large magnitude. 
The error was found to come from an improper use 
of the interpolation information when computing pre- 
vious time level data for newly created fringe points. 
This was fixed by changing the time-integration for 
such points to first-order. When updating the bound- 
ary points in the overset grids is done properly, the 
computed results are in fairly good agreement with the 
experimental data. However, the overlapped grid cal- 
culations still tend to generate erroneous oscillations 
in time. The composite grid computations also show 
some of this type of error, but to a much smaller ex- 
tent. Further investigation into the source of these 
errors is on going. 
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